Geometrical properties of the potential energy of the soft-sphere binary mixture 
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We report a detailed study of the stationary points (zero-force points) of the potential energy 
surface (PES) of a model structural glassformer. We compare stationary points found with two 
different algorithms (eigenvector following and square gradient minimization), and show that the 
mapping between instantaneous configuration and stationary points defined by those algorithms is 
as different as to strongly influence the instability index K vs. temperature plot, which relevance 
in analyzing the liquid dynamics is thus questioned. On the other hand, the plot of K vs. energy 
is much less sensitive to the algorithm employed, showing that the energy is the good variable to 
discuss geometric properties of the PES. We find new evidence of a geometric transition between 
a minima-dominated phase and a saddle-point-dominated one. We analyze the distances between 
instantaneous configurations and stationary points, and find that above the glass transition, the 
system is closer to saddle points than to minima. 



I. INTRODUCTION 

Interest in the geometric properties of the potential 
energy surface (PES) of liquids as a mean to under- 
stand their dynamics and thermodynamics dates back 
to the work of Goldstein [1] and Stillinger and Weber [2]. 
These works showed that to analyze the low temperature 
dynamics of supercooled liquids and glasses it is useful 
to map instantaneous configurations (ICs) to the (local) 
minimum of the PES found directly downhill, called in- 
herent structure (IS). This mapping allows to disregard 
fast vibrations, focusing on slow, activated, structural 
relaxations. But if one aims to describe the dynamic 
crossover taking place around the the Mode Coupling 
critical temperature (Imc) [3], the IS mapping is not 
useful because barriers between minima are no longer 
relevant at high temperatures and the two timescales 
cease to be well separated. Within the PES approach, 
two approaches have been proposed. One is to consider 
whole superstructures of minimia (called metabasins [4- 
6]). Another is to focus on stationary points with some 
unstable directions: saddle points (SPs) [7]. 

The latter approach was motivated by results obtained 
on the p-spin model (a mean field glass model). In this 
system, a threshold energy exists which separates a low- 
energy, minima-dominated region, from a high-energy, 
saddle dominated one [8, 9]. The higher the energy, the 
larger the number K of unstable directions of the typ- 
ical SP (this number is called the order, or instability 
index of the SP, and is equal to the number of negative 
eigenvalues of the Hessian matrix evaluated at the SP). 
In this model the average index as a function of the en- 
ergy K{E) can be computed [8]. Furthermore, it can be 
verified directly that at high temperatures the stationary 
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point closer to the typical IC is a SP with extensive K, 
while it is a minimum (K = 0) at low temperatures [10]. 
The dynamic arrest observed at the dynamic transition is 
related to the fact that the system starts getting trapped, 
or nearly trapped [11, 12]. A similar scenario was then 
proposed [7] for structural glasses: the glass transforma- 
tion is the consequence of a geometrical transition. In the 
saddle-dominated (high energy) phase, the system can re- 
lax either by jumping an energy barrier or by finding an 
unstable direction. In the minima-dominated phase (low 
energy), the second mechanism is no longer available; as 
a consequence relaxation times soar. It was shown later 
that the two-step relaxation observed in the supercooled 
liquid (and described by Mode Coupling theory (MCT)) 
can be qualitatively understood as relaxation in the vicin- 
ity of a SP [13]. 

This scenario has been explored in several numerical 
studies of model liquids. Some of these works obtained 
estimations of the K(E) curve [14-16], which have been 
found compatible with the existence of a geometric tran- 
sition (further evidence for the transition has been found 
in the context of high frequency vibrations [17]). Other 
works have studied instead A" as a function of the temper- 
ature T [18-22], showing that K decreases dramatically 
on approaching Tmc- I n early studies the view was held 
that a sharp transition can be observed as a function 
of temperature, with K = for T < Tmc and K > 
for T > Tmc, but further work has shown [5, 21, 23] 
that although K decreases very fast (most likely with an 
Arrenhius law [5, 24]), it is still nonzero for T < Tmc- 
This has prompted criticism of the saddle-minima transi- 
tion point of view (see e.g. ref. 24), although a geometric 
transition, controled by the energy, is compatible with a 
smooth K(T) curve (see sec. IV). 

But there is another issue to be discussed when con- 
sidering K(T) curves. Since the system is never precisely 
at a SP, to define a K(T) curve one needs to introduce 
a mapping between ICs and SPs (in a sense defining a 
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"basin of attraction" of a SP, and generalizing the IS con- 
cept). Additionally, if one wants to somehow interpret 
dynamic behavior from such curve, the mapping should 
preserve at least some dynamic information. In analyti- 
cal studies (e.g. [10, 25, 26]) ICs are (reasonably) mapped 
to the nearest SP (using the Euclidean distance or some 
overlap function). In contrast, in numerical work ICs 
are mapped to a SP through the algorithm used to find 
the latter, thus in principle introducing a dependence on 
the details of the procedure used to find SPs [25, 26] 
and raising the question of the dynamical relevance of 
the mapping. This is perhaps more worrying given that 
at least one popular procedure fails rather often, leaving 
some configurations unmapped (see ref. 21 and discus- 
sion below). Also puzzling are some results [21, 22] that 
seem to indicate that in the typical distance from an in- 
stantaneous configuration to a saddle or to a minimum 
is the same, at variance with the mean field situation. 

In this paper we address the issue of the mapping be- 
tween ICs and SPs, and analyze of distances between ICs 
and SPs and minima in more detail than has been pre- 
viously done. Our results show that the K vs. T plots 
are algorithm-dependent, and that, at least in the soft- 
sphere model we consider, ICs at high temperature are 
closer to SPs with K > 0. The large number of SPs col- 
lected allows a new analysis which provides new evidence 
for the existence of a geometrical transition. 



II. MODEL AND ALGORITHMS 



We have considered the soft sphere binary mixture 
[27, 28], which consists of 50% of particles of type A and 
50% of type B, interacting with a pair potential Vij(r) = 
(<Ji + <Tj) 12 /r 12 . The radii <7j are fixed by the conditions 
a B /(TA = 1.2 and (2a A ) 3 + 2{a A + a B f + (2a B ) 3 = 4. 
We have used a system of N = 70 particles at unit den- 
sity and a smooth (cubic polynomial) long-range cut-off 
at -\/3 as in ref. 16. We have used swap Monte Carlo 
[29] to equilibrate the system at temperatures T = 1, 
0.683, 0.482, 0.350, 0.260, 0.220. For this system T MC is 
about 0.24 [16, 28]. At each temperature, 40000 equili- 
brated configurations were saved and used as starting 
point for minima and saddle point searches. Minima 
were obtained with Nocedal and Liu's LBFGS algorithm 
[30], which code can be obtained from the internet [31]. 
For SP searches, two different algorithms were employed: 
square-gradient minimization (SGM) and eigenvector fol- 
lowing (EF) (described below), to compare two different 
IC-SP mappings. In all, about 3.2 • 10 5 SPs were ob- 
tained. 



A. Square gradient minimization 

One way of finding SPs is minimizing the squared mod- 
ulus of the gradient, 

N 3 / q t/ \ 2 
i— 1 a— 1 x ■ / 

Since 4> is nonnegative, at the absolute minima <j> = 0, 
which implies \7V = (a saddle point). This method is 
relatively easy to implement, since good numerical mini- 
mization algorithms are publicly available (we have used 
LBFGS [30, 31]). The biggest drawback is that mini- 
mization can (and does rather often) converge to a local 
minimum, which is neither a saddle point, nor close to 
one in any reasonable sense [32] . 

B. Eigenvector following 

This method has been specifically designed to find sta- 
tionary points of the potential energy. Based on an orig- 
inal proposal by Cerjan and Miller [33], it has been sub- 
stantially improved by others (see ref. 34 and references 
therein). The problem originally considered [33] was to 
find a saddle point of index 1 starting from a local min- 
imum of V. The idea was to consider the function on 
a small sphere around the minimum. Using Lagrange 
multipliers and a quadratic approximation, one looks for 
local minima of the function constrained to the sphere. 
Close enough to the initial point (minimum), there is 
one local minimum of the constrained function that has 
higher energy: this is a point along the path that leads 
to the sought saddle point, and is taken as the starting 
point of the next iteration. Close to the saddle this cri- 
terion no longer applies, so a Newton-Raphson step is 
taken. 

We have used our own implementation of the eigenvec- 
tor following method as described by Wales and cowork- 
ers [22, 35-37]. At each iteration a step Ax is proposed, 
which in the base that (locally) diagonalizes the Hessian 
is [35, 36] 

Ax, = S, 2g " , (2) 

im(i + v i+4 ^/^) 

where ft.„ are the eigenvalues of the Hessian and g^ are 
the components of the gradient in the diagonal base (Aa;^ 
is set to for the directions where = 0; i.e. uniform 
displacements). The sign = ±1 is chosen as explained 
below. Note that as — ► 0, 

= -^ + 0{g% g, (3) 

where the first term is the Newton-Raphson step. A set 
of trust radii {8^} is mantained (one for each direction) 
[37]. The proposed step is rescaled so that |Aa; M | < <5 M for 
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all fx, and then the position is updated. Initially, the 5^ 
arc set to 0.2, and at each step are increased (decreased) 
by a factor 1.2 according to whether the quantity r — 
(h e — h^/hfj, is less (greater) than 1. h e is an estimation 
of the eigenvalue, h e = (g^ — g'^/Ax'^, where the prime 
means the quantity evaluated at the previous iteration 
[37]. 

If = 1, the step increases the energy along the di- 
rection ll, thus the algorithm converges to a maximum 
along this direction. Conversely if = — 1 the algrithm 
converges to a minimum along ll. Since a saddle point of 
order if is a maximum along K directions and a mini- 
mum along 3iV — K directions, in princple the algorithm 
may be made to converge to a saddle of the desired in- 
dex by setting S M = — 1 for 1 < ll < K, and 5^ = 1 for 
ll > K. In this work we do not want to fix the index from 
the start of the search, so for each starting configurations 
we run 20 steps with = — sgn/i M and only then fix the 
index to whatever value it has reached after the first 20 
steps [22]. 
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FIG. 1: Root-mean-square gradient vs. energy of the configu- 
ration obtained with SGM or EF starting from 10000 equilib- 
rium configurations at T — 0.350. Inset: success rate (frac- 
tion of initial configurations that converge to saddle points) 
vs. temperature (SGM: circles, EF: squares). 



The distance we report is the Euclidean distance 



d= j£2{xi, a -yi, a ) 2 , (4) 

minimized over the symmetry operations of the system 
(i.e. translations, the 48 discrete symmetries of the simple 
cubic lattice, and particle permutations). Minimization 
over translations is done applying Brent's method [38, 
§10.2] to a distance minimized over the discrete symme- 
tries and permutations. This in turn is found by exhaus- 
tive exploration of the discrete symmetries and using the 
Hungarian algorithm [39] (as implemented by B. Gerkey 
[40]) to minimize over permutations. 

III. COMPARISON OF SGM AND EF 
MAPPINGS 

A. Basins of attraction and success rate 

In Fig. 1 where we show the rms gradient g = a/ <j)/N of 
the configurations obtained after running SGM and EF 
on the same set of ICs. EF produces tightly converged 
saddle points (g < 10~ n ), while the configurations found 
by SGM have rms gradients that cluster around about 
10~ 6 and 10~ 2 . To decide whether these configurations 
are saddle points and compute a success rate (Fig. 1, 
inset) we have used the SGM configurations as starting 
points for EF searches and computed the distance be- 
tween the SGM and corresponding EF converged config- 
urations. In some cases, after a few (3-5) steps, EF found 
a saddle very close to the SGM connfiguration (distances 
of the order of 10~ 5 to 10~ 2 ), while in other cases the 



distance was 0(1) or larger (and the number of steps 
increased). The first case clearly corresponds to an abso- 
lute minimum of <j>, and the second to a local minimum. 
The criterion we used to accept the SGM configuration 
as a true saddle was to require that the distance between 
it and the corresponding EF saddle be less than 0.01. 
This is approximately equivalent to requiring g < 10~ 4 
for the SGM configuration. 

The high failure rate of the SGM algorithm makes it 
unsuitable to define basins of attraction of saddle points 
[21, 22]. This failure rate is not due to problems with 
the minimization algorithm, but to the high number of 
local minima of the function <f> [14, 16, 21]. On the other 
hand the EF algorithm in principle will always converge 
to a saddle point (eventual failures being due to numer- 
ical problems or implementation details). However, it 
should be remembered that although basins of attraction 
for saddle points can be defined using EF, these are not 
necessarily a reasonable generalization of ISs. It is well 
known that iterative non-linear algorithms can lead to 
multiply connected or even fractal basins (a case of frac- 
tal basins is the Newton-Raphson algorithm applied to 
finding the roots of the polynomial z 3 — 1, see e.g. ref. 38, 
§9.4). In the case of EF, a detailed study on a 3-atom 
cluster [35] has shown that the basins, though not fractal, 
are still complex and multiply connected. Their relevance 
to liquid dynamics is thus not to be taken for granted. 
We have not performed such detailed analysis here, but 
in Fig. 2 we plot the instantaneous energy along a short 
(1000 steps) molecular dynamics (MD) run, along with 
the energy of the EF saddle points found starting from 
each IC. All these ICs map to a single IS (minimum). 
The strong energy fluctuations of the saddles found in 
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0.4 0.6 

MD time 

FIG. 2: Instantaneous energy along a 1000-step MD run at 
T — 0.260 (dashed-dotted line) and energy of the saddle 
points found by EF using the MD configurations as starting 
points (full line) . All instantaneous configurations correspond 
to a single IS. 



this way indicate that the usefulness of the basins of at- 
traction so defined is probably rather limited in under- 
standing the liquid dynamics. 



B. Saddle index curves 

Let us first consider the saddle index vs. temperature 
curve. In Fig. 3 we plot separately the average value 
of K for saddles obtained with EF and SGM. We also 
plot K(T) evaluated for ICs (the eigenvectors of the 
Hessian evaluated at an IC are usually called instanta- 
neous normal modes (INM)). The fact that the curves 
are algorithm-dependent prevents one from drawing any 
dynamical conclusion from them, unless there is some 
reason to expect that the mapping between ICs and SPs 
preserves some dynamical information. In particular, the 
critical temperature To where K(Tq) = (which might 
or might not be greater than zero) will be algorithm- 
dependent and thus not of much significance without fur- 
ther evidence of the dynamical relevance of the algorithm 
chosen to compute K(T). 

On the other hand, if we consider K as a funcion of 
the energy of the saddle point (Fig. 3, right), the curves 
produced from the SPs collected with EF and SGM are 
essentially coincident in the region of energies where both 
algorithms find a significant number of saddles (see in- 
set of Fig. 4). The corresponding curve for INMs (not 
shown) is very close to those corresponding to the SPs, 
in contrast to what is found in Lennard- Jones [14]. In 
this purely geometrical plot, the problem of the IC-SP 
mapping is avoided, and issues such as the existence of a 
geometrical transition can be meaningfully discussed. 
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FIG. 3: Average instability index vs. temperature (left) and 
energy (right) for INM, SGM saddles, and EF saddles. Each 
point is an average over SPs obtained from 40000 ICs. 



4 -i 



O 



1 - 




1.6 2.0 2.4 2.8 3.2 

Energy / N 



10 



K 



FIG. 4: Logarithm of the number of SPs found vs. instability 
index, for the energy band 1.9 < e < 2.0. Inset: Logarithm 
of the number of SPs vs. energy. 



Of course, this does not mean one should not worry 
about possible biases introduced by the algorithms. In- 
deed, if one considers the SPs within a given energy band, 
the distribution of the saddle index is slightly different, 
with EF tending to be slightly narrower (Fig. 4 shows a 
representative energy band). It is also clear that SGM 
tends to find SPs with lower energy (inset of Fig. 4). 
However, we find that the maximum of the logiV samp vs. 
K curves are the same for both algorithms at all energy 
bands. The significance of this maximum can be appre- 
ciated from the considerations of sec. IV. 
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FIG. 5: Average distance from instantaneous configurations 
vs. temperature for minima, EF saddles, and SGM saddles. 
Error bars are estimates of sample standard deviation (not 
errors on the average themselves). 



C. Distances 

Finally we consider the distances between ICs and the 
stationary points (IS, EF saddles and SGM saddles) ob- 
tained starting from the given IC. The distances reported 
here are those obtained after minimizing over the sym- 
metries of the hamiltonian, as discussed in sec. II C. We 
have found that this distance coincides most of the time 
with the distance obtained without applying minimiza- 
tions when one considers an IC and the stationary point 
obtained from it, so that the averages we report are not 
significantly different from those obtained without min- 
imizing. Minimization is however important when com- 
puting distances between a configuration and a station- 
ary point obtained from a different IC, as we do below. 

The average distances as a function of temperature are 
plotted in Fig. 5. Two things are to be noted: first, SGM 
saddles are always closer than EF saddles; second, ISs are 
farther than SGM saddles at high temperatures but start 
to be found closer as temperature is lowered. The first 
fact points to the influence of the algorithm in defining an 
IC-SP mapping, stressing the problems of interpretation 
of K(T) curves. The second provides direct evidence that 
at high temperatures there is a saddle point to be found 
closer to the typical IC than the corresponding IS. 

To investigate this matter more closely, we look in de- 
tail at a short MD trajectory at T = 0.26, slightly above 
Tmc (Fig. 6). For all ICs in this run, apart from comput- 
ing the corresponding IS and SGM saddle as usual, we 
have searched for the nearest SP in the pool of all SPs 
found at the corresponding temperature. We find again 
that the system is mostly close to a saddle point of order 
K > 1, as can be also seen in the inset. This result is 
natural within the geometrical transition scenario, but 



FIG. 6: Distance from instantaneous configurations along a 
short MD run at T = 0.26 to the corresponding IS, SGM 
saddle and nearest SP found. Inset: index of nearest SP. 

this is, to our knowledge, the first direct observation of 
this fact in a liquid. We furthermore find that the SGM 
saddle is not the closest saddle point. Of course, our pro- 
cedure does not guarantee to give the closest saddle to a 
given IC, but we do find SPs closer than the SGM saddle. 
We must note that Wales and Doyle [22], in an analysis 
of a Lennard- Jones binary mixture, did not find differ- 
ences in the mean distances from ICs to SPs or ISs. It 
would be interesting to analyze that system at the single 
configuration level (in the spirit of Fig. 6). 

IV. SADDLE-MINIMA TRANSITION 

Consider the number Af(K, E) of saddle points of or- 
der K and energy between E and E + dE. For short 
range interactions, one expects to be able to divide the 
system into effectively independent subsystems, so that 
Af(K, E) should be exponential in the size of the system 
[22, 41]. Then log N is extensive in the thermodynamic 
limit, and the complexity E(fc,e) = (1/N) \ogM{K, E) 
is an intensive quantity. The (intensive) average saddle 
index can then be written (e = E/N) 

(k(e)) = ±£°^exp[NX{K/N,e)]dK (5) 
N f°° 

= — fcexp[7V£(fc,e)]dfc, (6) 
Z Jo 

where 

Z= / exp[iVE(fc, e)] dk. (7) 
Jo 

Using the saddle point method one gets 

(k(e)) = k(e) + 0(l/N), (8) 
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FIG. 7: Logarithm of number of saddle points. 



where k(e) is the point where E(fc, e) attains a maximum 
(with e fixed), i.e. the solution of dT,(k,e)/dk = 0. The 
saddle-minima transition should be understood as hap- 
pening at the threshold energy eth, defined as the maxi- 
mum energy for which k(e) = 0. In the thermodynamic 
limit this implies (k(e < e t h)) = 0. In finite systems the 
average will remain positive, but the (k(e)) curve will 
show a fast crossover, remnant of the sharp N — > oo tran- 
sition, just as in other thermodynamic phase transitions. 
Clearly, the transition does not mean that there are no 
saddles for e < e t h (quite the contrary, there is an expo- 
nential number of them), but that they are subdominant 
respect to minima: Af(K > 0)/A/"(0) — > exponentially 
with N. 

The control parameter of the (geometric) saddle- 
minima transition is the energy [7, 14--17]. If one consid- 
ers K{T) (assuming one defines a mapping free from the 
problems discussed above), one will very likely find K > 
below Tmc because relaxation processes, though slow, 
will eventually sample saddle points (transition states). 
A smooth K(T) curve is thus compatible with a geomet- 
ric transition. Within the landscape point of view, one 
can understand the sharp dynamic crossover happening 
in fragile liquids around Tmc as the signature of a geo- 
metric transition [16, 17]. 

Given the large number of SPs obtained for this work 
(about 3.2- 10 5 ), we can try to obtain a rough estimate of 
the qualitative behavior of £(fc,e). The SPs have been 
classified into energy bands of width 0.1, and for each 
band a histogram in K was constructed. The logarithm 
of the histogram heights gives an estimate of the shape 
of the actual £ and is shown in Fig. 7. One can clearly 



see a maximum that goes to k = for low values of 
the energy. From the present data one obtains eth = 
1.77 ± 0.01, which is compatible with the values found in 
refs. 16 and 17. 

V. CONCLUSIONS 

We have shown that in numerical studies of the PES of 
liquids, the algorithm chosen to associate instantaneous 
configurations and saddle points can have a significant in- 
fluence in the analysis of quantities like the saddle index 
vs. temperature curve. Of the two algorithms used, we 
have found that the saddles found with SGM are closer to 
the instantaneous configurations than those found with 
EF. Since no such difference was found in ref. 22, the 
present results may reflect a property of the soft-sphere 
model, or of the present implementation of the EF algo- 
rithm. In any case, the point is that curves such as K(T) 
are not meaningful unless the validity of the IC-SP map- 
ping with respect to dynamic properties is established. It 
seems that neither of these algorithms is useful to define 
a partition of phase pase into basins of attraction of a SP 
(in the case of SGM it is rigorously impossible [21]). 

The natural variable to analyze geometrical proper- 
ties of the PES is the energy. We have shown that this 
choice of variable largely avoids the issue of the IC-SP 
mapping (in particular the K{E) plot is mostly indepen- 
dent of the mapping), though the algorithms introduce 
some detectable bias in the sampling. We have produced 
an estimate of the shape of the saddle complexity that 
provides new evidence for the existence of a geometric 
transition in the soft-sphere model. 

Finally, our analysis of distances has shown that in 
this model above Tmc the system is closer to saddle 
points than to inherent structures, as has been shown 
to be the case in some mean- field models (p-spin [10], 
fc-trigonometric [25]). We have also found that there are 
saddles closer than the SGM saddle (this has been ex- 
plored analytically in the /c-trigonometric model, where 
the SGM saddle was found to be the closest saddle [25], 
and the mean-field (j) 4 model, where it is not [26]). 
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